
function[gmm,W,g]=gmm_fun_y_c_avg_eta_H(candidate_params, candidate_indices, parameters,Yt, vdjQ, t, side, retailer, step)

global  xij_for_investor eta Dmax Imax homedealer

%parmaeters
parameters(candidate_indices)=candidate_params; 

%sort Yt on id and then dealer 
if homedealer==0
    Yt = sortrows(Yt, [1,16,17]); 
else
    Yt = sortrows(Yt, [1,16,26]); 
end 

%number of trades    
Ntrades  = size(Yt,1);

%index for dealers
if homedealer==0 %benchmark model with no homedealers
    idxj        = find(Yt(:,17)==1); %index with first entry of each dealer 
    dealer_ids  = Yt(idxj,16);
    idxd        = find(ismember(Dmax,dealer_ids));
else %extension with  homedealers
    idxj        = find(Yt(:,17)==1);
    dealer_ids  = Yt(idxj,16);
    idxd        = find(ismember(Dmax,dealer_ids));
    
    idxjH        = find(Yt(:,26)==1 & Yt(:,25)==1);
    dealer_idsH  = Yt(idxjH,16);
    idxdH        = find(ismember(Dmax,dealer_idsH));
    
    idxjB        = find(Yt(:,26)==1 & Yt(:,25)==0);    
    dealer_idsB  = Yt(idxjB,16);
    idxdB        = find(ismember(Dmax,dealer_idsB));
end 

%label parameters
sigma   = Yt(1,15); 

if homedealer==0
    xij     = Yt(idxj,14);

else
    Jt = size(idxj,1);
    
    xijH        = nanmedian(Yt(idxj,24))*ones(Imax,1); 
    xijH(idxd) = Yt(idxj,24); 

    xijB        = nanmedian(Yt(idxj,23))*ones(Imax,1); 
    xijB(idxd) = Yt(idxj,23); 
    
    %TO RUN THE EXTENDED MODEL, SETTING HOMEDEALER BENEFIT TO 0, USE THE BELOW
    %xijH        = nanmedian(Yt(idxj,23))*ones(Imax,1); 
    %xijH(idxd) = Yt(idxj,23); 

    %xijB        = nanmedian(Yt(idxj,23))*ones(Imax,1); 
    %xijB(idxd) = Yt(idxj,23); 

    %set the quality of the dealer who is inactive to NaN (if Jt<Imax)
    xijB(setdiff(1:Imax, idxd))=NaN;
    xijH(setdiff(1:Imax, idxd))=NaN;
    
    %create a xij matrix with xij's from each dealer j's percpective
    xij = NaN*ones(Imax,Imax); 
        %Note: matrix with the xij's from dealer j's perspective
                   %column j are the xij's from dealer j's perspective,
                   %with (j,j) as her home dealer quality and all other
                   %entries as non-homedealer qualities
    for j=1:Imax
        xij(:,j) = xijB;
        xij(j,j) = xijH(j); %quality with base qualities of all dealers buy j and homedealer quality for j
    end
    
    %set quality of the dealer who is inactive to NaN (if Jt<Imax)
    xij(:, setdiff(1:Imax, idxd))=NaN;
end   

if xij_for_investor==1         %Benchmark model - investor gets 0 surplus
    xid     = 0;                      %set to 0 if psi does not include xid
    xij_shocks = Yt(:,14) ;    %set to 0 if price excludes xij 
    
elseif  xij_for_investor==0 
    %Note: This is a robusntess check where we leave xij to the investor as a surplus when eta=0 and homedealer==0
    
    if homedealer==0
        xid        = Yt(idxj,14);  
        xij_shocks =  0; 
        
    elseif homedealer==1 && eta==0
        xid = Yt(idxj,24) -Yt(idxj,23);        %homedealer benefit
        
        %TO RUN THE EXTENDED MODEL, SETTING HOMEDEALER BENEFIT TO 0, USE THE BELOW
        %xid = (Yt(idxj,24) -Yt(idxj,23))*0;  %JUST FOR CHECKING
        
        xij_shocks = Yt(:,23);                  %base quality

    end 
end 

q       = Yt(idxj,18);

if retailer==0
    munu    = parameters(1);
    sigmanu = parameters(2);
    c       = parameters(3);
elseif retailer==1
    munu    = parameters(1);
    sigmanu = parameters(2);
end 

%distribution functions
f = @(x) normpdf(x,munu,sigmanu);
F = @(x) normcdf(x,munu,sigmanu);


%Back out shocks nu    
theta_long   = Yt(:,5); %midbloomberg price
price            = Yt(:,8); %OTC price 

if  eta==0
    if side==1
        shocks = (price - theta_long) + xij_shocks;
    elseif side==-1
        shocks = (price - theta_long) - xij_shocks;
    end

elseif  eta~=0 && xid==0 && retailer==1 && homedealer==0 
    %Note: this is only done for the benchmark model, where the investor gets 0 surplus OTC and w/o homedealer
   
    vdjQ_long = zeros(size(Yt,1),1);
    for t=1:size(Yt,1)  %Note: careful this only works if Yt is sorted!
        idx0            = find(Yt(t,16)==dealer_ids); 
        vdjQ_long(t)    = vdjQ(dealer_ids(idx0));
    end 

    if side==1 %price = eta*value + (1-eta)(theta + nu - xij_shocks) -> nu = (price - eta*value)/(1-eta) - theta + xij_shock
        shocks = (price -eta*vdjQ_long)./(1-eta)- theta_long + xij_shocks;
    elseif side==-1
        shocks = (price -eta*vdjQ_long)./(1-eta) - theta_long - xij_shocks;
    end
end 
    


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INSTITUIONAL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if retailer==0
    
%use mean theta per dealer (to have a single theta per day)
thetaj = zeros(Imax,1);
for i=1:Imax
    thetaj(i)  = mean(Yt(Yt(:,16)==Dmax(i),5));
end 

%If eta!=0, compute the value of the dealer for the given parameters for benchmark model
if eta~=0 && xid==0 && retailer==0 && homedealer==0
     %Note: this is only done for the benchmark model, where the investor gets 0 surplus OTC and w/o homedealer
      estimate=1;
     [~,~,~,~,~,~,~,vdj,~,~,~,~,~,~,~,~,~,~,~,~,~,~,~] ...
         = model_prediction_avg_eta(q, Yt, parameters, side, estimate, 1,[],idxj, 0); 
end 

%initialize so that all outputes have the same number of dealers
psi        =  NaN*ones(Imax,1);

if side==1
    if homedealer==0
        psi(idxd,:) =  sigma*log(sum(exp(1/sigma*(xij+q))))-xid;
    else
        for j=1:Jt
            if xij_for_investor==1
                    xidj=0;
                else 
                    xidj=xid(j);
            end
            psi(idxd(j),:) =  sigma*log(sum(exp(1/sigma*(xij(idxd,idxd(j))+q))))-xidj;
        end
    end 
    
elseif side==-1
    if homedealer==0
        psi(idxd,:) =  -sigma*log(sum(exp(1/sigma*(xij-q))))+xid;
    else  
        for j=1:Jt
             if xij_for_investor==1
                    xidj=0;
                else 
                    xidj=xid(j);
             end 
            psi(idxd(j),:) =  -sigma*log(sum(exp(1/sigma*(xij(idxd,idxd(j))-q))))+xidj;
        end
    end 
end 

%into long shape 
psi_long   = zeros(size(Yt,1),1);
theta_long = zeros(size(Yt,1),1);
vdjQ       = zeros(size(Yt,1),1);

for t=1:size(Yt,1)  %Note: careful this only works if Yt is sorted!
    idx0            = find(Yt(t,16)==dealer_ids); 
    psi_long(t)     = psi(dealer_ids(idx0));
    theta_long(t)   = thetaj(dealer_ids(idx0));   
    
    if eta~=0 && xid==0  && homedealer==0
    vdjQ(t)         = vdj(dealer_ids(idx0));
    end 
end 

if eta~=0 && xid==0  && homedealer==0
    if side==1 %price = eta*value + (1-eta)(theta + nu - xij_shocks) -> nu = (price - eta*value)/(1-eta) - theta + xij_shock
    shocks = (price -eta*vdjQ)/(1-eta)- theta_long + xij_shocks;
    elseif side==-1
    shocks = (price -eta*vdjQ)/(1-eta) - theta_long - xij_shocks;
    end
end 
    
  
if side==1
      
        b=Inf;
        a      = (psi_long + c - eta *(vdjQ + xij_shocks))./(1-eta) -theta_long;
           
        Mills   = f(a)./(1 - F(a));  
        ExpfunT = munu + sigmanu.*Mills ;  
        VarfunT = (sigmanu.^2 .*(1 + (a .*f(a))./(1 - F(a)) - (Mills).^2));

        massOTC  = F(b) - F(a);
        massE    = F(a);
        
elseif side==-1 
        
        b = (psi_long - c + eta *(vdjQ + xij_shocks))./(1-eta) -theta_long ; 
        a= -Inf;
        
        Mills   = - f(b)./F(b);  
        ExpfunT = munu + sigmanu.*Mills ; 
        VarfunT = (sigmanu.^2 .*(1 + (-b .*f(b))./(F(b)) - (Mills).^2));

        massOTC  = F(b) - F(a);
        massE    = 1-F(b);
end   
    
        probOTC  = massOTC./(massOTC+massE); %probability to choose OTC conditional con trading      
        shareOTC = size(Yt(Yt(:,9)==1,:),1)./Ntrades;
       

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% RETAILER
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elseif retailer==1
    
    if side==1
        
        ExpfunT = munu;
        VarfunT = sigmanu.^2;
        
    elseif side==-1 
        
        ExpfunT = munu;
        VarfunT = sigmanu.^2;
    end  
    
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% MOMENTS AND GMM FUNCTION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

moment1 = (ExpfunT - shocks);
moment2 = (VarfunT -(shocks -ExpfunT).^2);

if retailer==0
    
    idxE = find(Yt(:,9)==2);
    moment1(idxE)=0;
    moment2(idxE)=0;
    
    moment3 = probOTC - shareOTC;
    moment = [moment1, moment2, moment3];
    
elseif retailer==1
    moment = [moment1, moment2];
end 

g =  1/Ntrades *sum(moment);

if step==1 
    W=eye(size(g,2));
elseif step==2 
    W=inv(1/(Ntrades) * (moment') * moment);
end 

gmm0 = g* W* g';

if isnan(gmm0)==0 
        gmm = gmm0;
   else 
        gmm = 10^100;
end 
    
  






